ESA18 TERN Workshop: Landscape Assessment (AusCover) - Remote Sensing

 

"Simple Change Detection in Raster using TERN web services"

     

This document contains the Landscape Assessment (AusCover) section of TERN’s workshop at the Ecological Society of Australia 2018 conference. It is based on Peter Scarth’s tutorial "Simple Change Detection in Raster using TERN web services and Google Earth Engine". Although it (more or less) follows Dr. Scarth’s tutorial steps, it has some variation in its content. Dr. Scarth’s tutorial uses a Jupyter Notebook with code in Python. In the workshop we will use R.

To analyse changes in Normalized Difference Vegetation Index (NDVI) over time, we will use datasets containing composited seasonal surface reflectance images (4 seasons per year). The datasets were created from full time series of Landsat TM/ETM+ imagery. The imagery has been composited over a season to produce new imagery that is representative of the period. In this process, techniques that reduce contamination by cloud and other problems were used.

In the final part of this section of the workshop we explore the causes of NDVI change in the area of most noticeable greening. To do so, a time series of the green fraction is plotted. The green fraction is one of the bands in the fractional cover product. This is another AusCover product hosted by TERN.

In this section of the workshop, we will cover the following aspects:

  1. INTRODUCTORY STEPS (getting acquainted with using raster files in R):
  1. NORMALIZED DIFFERENCE VEGETATION INDEX (NDVI) CHANGE ANALYSIS:

For conciseness the messages returned by R are not displayed (i.e. only the results).

The images and plots in this report can appear small at times. However, they look fine on a computer monitor when the code is run in a computer.

INTRODUCTORY STEPS (getting familiar with the use of raster files in R)

Preparation: Getting ready

We start by: loading the required libraries, setting up a working directory, and cleaning up memory.

# Load  Libraries
# ===============

library(dplyr)
library(reshape2)
library(stringr)

library(RColorBrewer)
library(ggplot2)
library(gridExtra)
library(RStoolbox)

library(sp)

library(rgdal)
library(raster)
library(rasterVis)

##library(RgoogleMaps)
library(ggmap)

#library(maps)
#library(mapdata)
#library(maptools)


# Set Data Path (TERN AusCover data repository)
# =============================================
data.path = "http://qld.auscover.org.au/public/data/landsat/surface_reflectance/aus"



# Optional Steps (remove comments to run them)
# ==============

# Setting up my Current Working Directory
# ---------------------------------------
#getwd()
#my.CWDir = "C:/Users/uqbblanc/Documents/TERN/CWDir"
#setwd(my.CWDir)
#getwd()

# Clean up Memory
# ---------------
#rm(list=ls())
#list.files()

Downloading, Reading and Exploring raster datasets.

Next, we download the file containing composited seasonal surface reflectance images for Australia. The file contains 6 layers corresponding to different satellite bands. Here, using a vrt file, we load each band individually. ‘VRT’ files contain descriptions of >= 1 datasets, usually in XML format, to produce a virtual GDAL(Geospatial Data Abstraction Library) dataset composed from other GDAL datasets (see Part 2 of this tutorial and references within for further description). Typing the name of the raster object will provide it a summary of its features (dimensions, resolution, extent, coordinates reference system,..).

getwd()
## [1] "C:/Users/uqbblanc/Documents/TERN/03-Trainning_Conferences/181125_ESA18/03-Workshop/CWDir_Final/Final/2-AusCover/InGitHub_190205/New"
# Data Path and Name
data.path = "http://qld.auscover.org.au/public/data/landsat/surface_reflectance/aus"
data.fn = "l8olre_aus_m201609201611_dbia2.vrt"

# Download the data file (if it doesn't work try method='wget' or functions in Library 'RCurl')
download.file(url=paste(data.path, data.fn, sep="/"), destfile='ASLM.vrt', method='auto')
list.files()
## [1] "ASLM.vrt"                           
## [2] "ESA18_Workshop_AusCover_P1_v7f.html"
## [3] "ESA18_Workshop_AusCover_P1_v7f.Rmd"
# Load the data file => We see that it has 6 bands, but only one is loaded.
ASLM.ras = raster("ASLM.vrt") # See it has 6 bands
ASLM.ras
## class       : RasterLayer 
## band        : 1  (of  6  bands)
## dimensions  : 135159, 141481, 19122430479  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -1944645, 2299785, -4910195, -855425  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\uqbblanc\Documents\TERN\03-Trainning_Conferences\181125_ESA18\03-Workshop\CWDir_Final\Final\2-AusCover\InGitHub_190205\New\ASLM.vrt 
## names       : ASLM 
## values      : -32768, 32767  (min, max)
# Load the remaining 5 Bands.
ASLM.Band2.ras = raster("ASLM.vrt", band=2)
ASLM.Band2.ras
## class       : RasterLayer 
## band        : 2  (of  6  bands)
## dimensions  : 135159, 141481, 19122430479  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -1944645, 2299785, -4910195, -855425  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\uqbblanc\Documents\TERN\03-Trainning_Conferences\181125_ESA18\03-Workshop\CWDir_Final\Final\2-AusCover\InGitHub_190205\New\ASLM.vrt 
## names       : ASLM 
## values      : -32768, 32767  (min, max)
ASLM.Band3.ras = raster("ASLM.vrt", band=3)
ASLM.Band3.ras
## class       : RasterLayer 
## band        : 3  (of  6  bands)
## dimensions  : 135159, 141481, 19122430479  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -1944645, 2299785, -4910195, -855425  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\uqbblanc\Documents\TERN\03-Trainning_Conferences\181125_ESA18\03-Workshop\CWDir_Final\Final\2-AusCover\InGitHub_190205\New\ASLM.vrt 
## names       : ASLM 
## values      : -32768, 32767  (min, max)
ASLM.Band4.ras = raster("ASLM.vrt", band=4)
ASLM.Band4.ras
## class       : RasterLayer 
## band        : 4  (of  6  bands)
## dimensions  : 135159, 141481, 19122430479  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -1944645, 2299785, -4910195, -855425  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\uqbblanc\Documents\TERN\03-Trainning_Conferences\181125_ESA18\03-Workshop\CWDir_Final\Final\2-AusCover\InGitHub_190205\New\ASLM.vrt 
## names       : ASLM 
## values      : -32768, 32767  (min, max)
ASLM.Band5.ras = raster("ASLM.vrt", band=5)
ASLM.Band5.ras
## class       : RasterLayer 
## band        : 5  (of  6  bands)
## dimensions  : 135159, 141481, 19122430479  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -1944645, 2299785, -4910195, -855425  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\uqbblanc\Documents\TERN\03-Trainning_Conferences\181125_ESA18\03-Workshop\CWDir_Final\Final\2-AusCover\InGitHub_190205\New\ASLM.vrt 
## names       : ASLM 
## values      : -32768, 32767  (min, max)
ASLM.Band6.ras = raster("ASLM.vrt", band=6)
ASLM.Band6.ras
## class       : RasterLayer 
## band        : 6  (of  6  bands)
## dimensions  : 135159, 141481, 19122430479  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : -1944645, 2299785, -4910195, -855425  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\uqbblanc\Documents\TERN\03-Trainning_Conferences\181125_ESA18\03-Workshop\CWDir_Final\Final\2-AusCover\InGitHub_190205\New\ASLM.vrt 
## names       : ASLM 
## values      : -32768, 32767  (min, max)

Subsetting rasters

The raster objects we have created contain seasonal surface reflectance images for the whole of Australia. They are very large objects, which makes their manipulation and processing slow. We now subset these objects to create smaller objects for our area and bands of interest.

# Create a new extent containing our area of interest
new.extent = extent(2030000, 2070000, -3160000, -3140000)
new.extent
## class       : Extent 
## xmin        : 2030000 
## xmax        : 2070000 
## ymin        : -3160000 
## ymax        : -3140000
# Crop (i.e. subset) the raster layers to our area of interest
ASLM.Band5SWIR.Bris.rL = crop(ASLM.Band5.ras, new.extent)
ASLM.Band5SWIR.Bris.rL
## class       : RasterLayer 
## dimensions  : 667, 1334, 889778  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 2029995, 2070015, -3159995, -3139985  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\uqbblanc\AppData\Local\Temp\RtmpmQk4lF\raster\r_tmp_2019-02-05_180608_14960_19577.grd 
## names       : ASLM 
## values      : -75, 7475  (min, max)
ASLM.Band4NIR.Bris.rL = crop(ASLM.Band4.ras, new.extent)
ASLM.Band4NIR.Bris.rL
## class       : RasterLayer 
## dimensions  : 667, 1334, 889778  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 2029995, 2070015, -3159995, -3139985  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\uqbblanc\AppData\Local\Temp\RtmpmQk4lF\raster\r_tmp_2019-02-05_180645_14960_01628.grd 
## names       : ASLM 
## values      : -47, 7554  (min, max)
ASLM.Band3Red.Bris.rL = crop(ASLM.Band3.ras, new.extent)
ASLM.Band3Red.Bris.rL
## class       : RasterLayer 
## dimensions  : 667, 1334, 889778  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 2029995, 2070015, -3159995, -3139985  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\uqbblanc\AppData\Local\Temp\RtmpmQk4lF\raster\r_tmp_2019-02-05_180712_14960_42801.grd 
## names       : ASLM 
## values      : 12, 7058  (min, max)
# Remove unneeded raster objects
rm(list = ls(pattern=".ras"))
ls()
## [1] "ASLM.Band3Red.Bris.rL"  "ASLM.Band4NIR.Bris.rL" 
## [3] "ASLM.Band5SWIR.Bris.rL" "data.fn"               
## [5] "data.path"              "new.extent"

Combining Raster Layers into a Multi-Layer Raster & Visualising the Composite Image

We now visualise the composite image of the 3 raster bands. To do so we will follow these steps:

  1. Remove (erroneous) negative values by giving them a missing value code (‘NA’)
  2. Create a multi-layered raster object from 3 single layered rasters. Specifically we combine rasterLayer objects into a rasterBrick object. See TERN DSDP "Using Raster Data in R" tutorial for more information on the raster objects classes in R.
  3. Plot the resulting rasterBrick using the plotRGB function in the raster package and the ggRGB function in the RStoolbox library with different parameters. The latter function uses the ggplot2 graphic library framework, so it is compatible qqplots2 features. That is, ggplot2 features can be added to ggRGB plots to enhance them.
# Replace values < 0 with NAs
# ===========================
 # Band 5: SWIR
 # ------------
ASLM.Band5SWIR.Bris.rL[ASLM.Band5SWIR.Bris.rL < 0] = NA
ASLM.Band5SWIR.Bris.rL
## class       : RasterLayer 
## dimensions  : 667, 1334, 889778  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 2029995, 2070015, -3159995, -3139985  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : ASLM 
## values      : 0, 7475  (min, max)
freq(is.na(ASLM.Band5SWIR.Bris.rL))
##      value  count
## [1,]     0 873072
## [2,]     1  16706
freq(is.na(ASLM.Band5SWIR.Bris.rL))/ncell(ASLM.Band5SWIR.Bris.rL)
##             value      count
## [1,] 0.000000e+00 0.98122453
## [2,] 1.123876e-06 0.01877547
 # Band 4: NIR
 # ------------
ASLM.Band4NIR.Bris.rL[ASLM.Band4NIR.Bris.rL < 0] = NA
ASLM.Band4NIR.Bris.rL
## class       : RasterLayer 
## dimensions  : 667, 1334, 889778  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 2029995, 2070015, -3159995, -3139985  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : ASLM 
## values      : 3, 7554  (min, max)
freq(is.na(ASLM.Band4NIR.Bris.rL))
##      value  count
## [1,]     0 873215
## [2,]     1  16563
freq(is.na(ASLM.Band4NIR.Bris.rL))/ncell(ASLM.Band4NIR.Bris.rL)
##             value      count
## [1,] 0.000000e+00 0.98138524
## [2,] 1.123876e-06 0.01861476
 # Band 3: Red
 # ------------
ASLM.Band3Red.Bris.rL[ASLM.Band3Red.Bris.rL < 0] = NA
ASLM.Band3Red.Bris.rL
## class       : RasterLayer 
## dimensions  : 667, 1334, 889778  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 2029995, 2070015, -3159995, -3139985  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : ASLM 
## values      : 12, 7058  (min, max)
freq(is.na(ASLM.Band3Red.Bris.rL))
##      value  count
## [1,]     0 873221
## [2,]     1  16557
freq(is.na(ASLM.Band3Red.Bris.rL))/ncell(ASLM.Band3Red.Bris.rL)
##             value      count
## [1,] 0.000000e+00 0.98139199
## [2,] 1.123876e-06 0.01860801
# Create a multi-layered raster object
# ====================================
ASLM.Bands3to5.Bris.rB = brick(ASLM.Band3Red.Bris.rL, ASLM.Band4NIR.Bris.rL, ASLM.Band5SWIR.Bris.rL)
names(ASLM.Bands3to5.Bris.rB) = c("Band3.Red", "Band4.NIR", "Band5.SWIR")
ASLM.Bands3to5.Bris.rB
## class       : RasterBrick 
## dimensions  : 667, 1334, 889778, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 2029995, 2070015, -3159995, -3139985  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : Band3.Red, Band4.NIR, Band5.SWIR 
## min values  :        12,         3,          0 
## max values  :      7058,      7554,       7475
# Plot the Composite Image of the 3 Raster Bands
# ==============================================

# Plot using function 'plotRGB' from 'raster' library
# ---------------------------------------------------
plotRGB(ASLM.Bands3to5.Bris.rB, r=1, g=2, b=3)  # 2nd best looking one

plotRGB(ASLM.Bands3to5.Bris.rB, r=3, g=2, b=1)  # Best looking one

# Plot using function ggRBG from package RStoolbox for more control and enhanced plots
# ------------------------------------------------------------------------------------
ggRGB(ASLM.Bands3to5.Bris.rB, r=3, g=2, b=1) + 
labs(title= "Landsat Surface Reflectance (Bands 5,4,3)", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=20, face="bold"), 
      axis.title = element_text(size=12, face="bold"), axis.text=element_text(size=9) )

# Stretch the values to increase the contrast of the image
# Histogram stretch is the nicest looking one (square-root and linear stretch not too bad; logarithmic stretch not good)
ggRGB(ASLM.Bands3to5.Bris.rB, r=1, g=2, b=3, stretch='hist') + 
labs(title= "Landsat Surface Reflectance (Bands 5,4,3)", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=20, face="bold"), 
      axis.title = element_text(size=12, face="bold"), axis.text=element_text(size=9) )

Re-project the Raster Object & Using a Map to find our bearings

We will use the functions in the package ggmap to query map servers for a map and plot the map. The function get_map in this package can query Google Maps, Naver Maps, OpenStreeMap, and Stamen Maps. Unfortunately, since July 2018 to download maps from Google Maps and Naver Maps an Application Program Interface (API) Key is required and there is an associated download costs. However, we can still freely download maps from Stamen.

The function ggmap from the ggmap package plots the raster objects produced by get_map. ggmap is also built ‘on top’ the ggplot2 graphic library framework, so can use this library features.

ggmap retrieves maps in geographic coordinates (latitude/longitude) on the World geodetic System of 1984 (WGS84) datum. Therefore, to compare our satellite images to these maps we first need to reprojected to this Coordinate Reference System (CRS).

CRS can be stored in many different formats. The most common ones include: PROJ.4 (the default output from many spatial data R packages; e.g. sp, rgdal, raster); EPSG codes, (a more concise format, found for example in the package sf); and Well Known Text (WKT) format. The EPSG format for geographic coordinates on the WGS84 datum is EPSG:4326. For more details on CRS and their formats, see TERN DSDP "Using Raster Data in R" tutorial.

We use grid.arrange from the package gridExtra place multiple ‘grobs’ (grid graphical objects, in our case ggplot2 derived graphics) in a page.

# View current CRS in PROJ.4 format
# ---------------------------------
crs(ASLM.Bands3to5.Bris.rB)  
## CRS arguments:
##  +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0
## +ellps=GRS80 +units=m +no_defs
# View EPSG codes as PROJ.4 formats
# ---------------------------------
 # The RasterBrick is currently in EPSG:3577. An adequate CRS for this location
CRS("+init=epsg:3577")
## CRS arguments:
##  +init=epsg:3577 +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0
## +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
## +units=m +no_defs
 # We want to re-projected to this CRS
CRS("+init=epsg:4326")
## CRS arguments:
##  +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
## +towgs84=0,0,0
# Reproject the RasterBrick & Create a Plot of the Re-projected RasterBrick
# -------------------------------------------------------------------------
ASLM.Bands3to5.Bris.rB.prjEPSG4326 = projectRaster(ASLM.Bands3to5.Bris.rB, crs="+init=epsg:4326")
#crs(ASLM.Bands3to5.Bris.rB.prjEPSG4326)
Bris.p1 = ggRGB(ASLM.Bands3to5.Bris.rB.prjEPSG4326, r=1, g=2, b=3, stretch='hist') + 
labs(title= "Landsat Surface Reflectance (Bands 5,4,3)", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
      axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) ) 

# Make a Bounding Box
# -------------------
 # Use extent of the re-projected RasterBrick
Bris.ext.GeoCoord = extent(ASLM.Bands3to5.Bris.rB.prjEPSG4326)
Bris.ext.GeoCoord
## class       : Extent 
## xmin        : 152.8747 
## xmax        : 153.3158 
## ymin        : -27.59609 
## ymax        : -27.35823
 # Create Bounding Box for Brisbane zone in Geographic Coordinates 
Bris.long = c(extent(Bris.ext.GeoCoord)[1:2])
Bris.lat = c(extent(Bris.ext.GeoCoord)[3:4])
bbox = make_bbox(Bris.long, Bris.lat, f=0.)
bbox
##      left    bottom     right       top 
## 152.87466 -27.59609 153.31583 -27.35823
# Get Map of Map from Stamen Maps & Create Plot of this Map
# ---------------------------------------------------------
#stamen.Bris.map = get_map(bbox, maptype="terrain", source="google") # Doesn't currently work
stamen.Bris.map = get_stamenmap(bbox, maptype="terrain")
Bris.p2 = ggmap(stamen.Bris.map) + labs(title= "Stamem Map", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
      axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) )

# Create a new smaller common extent to subset the image and map
# --------------------------------------------------------------
# The satelite images are oblique to the map. Focus on the shared extent. 
common.extent = extent(152.92, 153.25, -27.535, -27.425)
common.extent
## class       : Extent 
## xmin        : 152.92 
## xmax        : 153.25 
## ymin        : -27.535 
## ymax        : -27.425
# Crop & Create a Plot of the RasterBrick
# ---------------------------------------
# Crop to Raster Brick to common extent
ASLM.Bands3to5.Bris.rB.prjEPSG4326.common = crop(ASLM.Bands3to5.Bris.rB.prjEPSG4326, common.extent)
ASLM.Bands3to5.Bris.rB.prjEPSG4326.common
## class       : RasterBrick 
## dimensions  : 417, 1089, 454113, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 0.000303, 0.000264  (x, y)
## extent      : 152.9201, 153.2501, -27.53511, -27.42502  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : Band3.Red, Band4.NIR, Band5.SWIR 
## min values  : 54.181408, 53.980226,   6.303069 
## max values  :  6819.446,  7119.855,   6393.382
# Create Raster Brick Plot
Bris.p3 = ggRGB(ASLM.Bands3to5.Bris.rB.prjEPSG4326.common, r=1, g=2, b=3, stretch='hist') + 
labs(title= "Landsat Surface Reflectance (Bands 5,4,3) - Centre", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
      axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) ) 
      
# Get Map with the common extent & Create Plot of this Map
# --------------------------------------------------------
stamen.Bris.map.common = get_stamenmap(common.extent[c(1,3,2,4)], maptype="terrain")
Bris.p4 = ggmap(stamen.Bris.map.common) + labs(title= "Stamem Map - Centre", x="Longitude", y="Latitude") +
theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
      axis.title = element_text(size=11, face="bold"), axis.text=element_text(size=9) )
  
# Show All 4 Plots
# ----------------
# `Grid.arrange (from Package `gridExtra): Set up a layout to place multiple ggplots on a page.
grid.arrange(Bris.p1, Bris.p2, Bris.p3, Bris.p4, nrow=2)